from eofs.xarray import Eof
import xarray as xr
import numpy as np
import datetime as dt

import matplotlib.pyplot as plt  #用来绘图的包
import cartopy.crs as ccrs       
import cartopy.feature as cfeature 
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
import cmaps

file = xr.open_dataset('./mydata/sst_month_mean.nc')
file

<xarray.Dataset>
Dimensions:  (lat: 89, lon: 180, time: 2019)
Coordinates:
  * lat      (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
  * lon      (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
  * time     (time) datetime64[ns] 1854-01-01 1854-02-01 ... 2022-03-01
Data variables:
    sst      (time, lat, lon) float32 ...

sst = file.sst.groupby(file.time.dt.year).mean().rename({'year':'time'})
sst = sst.sel(time=slice('1979','2021'))
sst

<xarray.DataArray 'sst' (time: 43, lat: 89, lon: 180)>
array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        ...,
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]],

       [[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
...
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]],

       [[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        ...,
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
  * lon      (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
  * time     (time) int64 1979 1980 1981 1982 1983 ... 2017 2018 2019 2020 2021

距平值

# eof
solver = Eof(sst)
eof = solver.eofs(neofs=3,eofscaling=2)
pc = solver.pcs(npcs=3,pcscaling=1)
print('eof:',eof.shape)
print('-------------------')
print('pc:',pc.shape)

eof: (3, 89, 180)
-------------------
pc: (43, 3)

# 方差贡献率
var=solver.varianceFraction(neigs=3)
var

<xarray.DataArray 'variance_fractions' (mode: 3)>
array([0.32112256, 0.16426052, 0.06989542], dtype=float32)
Coordinates:
  * mode     (mode) int32 0 1 2
Attributes:
    long_name:  variance_fractions

year = np.arange(1979,2022)

def mapart(ax): #定义一个绘制地图地图的函数
    '''
    添加地图元素
    '''
    projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影，并指定中心经度为180°
    ax.coastlines(color='k', lw=0.5)                #添加海岸线
    ax.add_feature(cfeature.LAND, facecolor='white') #添加陆地
    # 设置地图范围，经度为（0，358），维度为（-88，88）
    ax.set_extent([0, 358, -88, 88], crs=ccrs.PlateCarree(central_longitude=180))
    # 设置经纬度标签
    ax.set_xticks([-120, -60, 0, 60, 120], crs=projection)
    ax.set_yticks([-80, -40, 0, 40, 80], crs=projection)
    lon_formatter = LongitudeFormatter(zero_direction_label=True) 
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)


def eof_contourf(EOFs, PCs, pers, year): 
    plt.close
    pers=(pers*100).values 
    fig = plt.figure(figsize=(18, 12), dpi=300)       
    projection = ccrs.PlateCarree(central_longitude=180) 
    ax1 = fig.add_subplot(3, 2, 1, projection=projection) 
    mapart(ax1) 
    p = ax1.contourf(file.lon,
                     file.lat,
                     EOFs[0], 
                     np.linspace(-1,
                                 1,
                                 41),
                     cmap=cmaps.BlueWhiteOrangeRed,
                     transform=ccrs.PlateCarree())  
    ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left') 

    ax2 = fig.add_subplot(3, 2, 2) 
    ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--') 
    b = ax2.bar(year, PCs[:, 0], color='r') 
    for bar, height in zip(b, PCs[:, 0]):
        if height < 0:
            bar.set(color='blue') 
    ax2.set_title('PC1' % (round(pers[0], 2)), loc='left') 

    ax3 = fig.add_subplot(3, 2, 3, projection=projection) 
    mapart(ax3)
    pp = ax3.contourf(file.lon,
                      file.lat,
                      EOFs[1],
                      np.linspace(-1,
                                  1,
                                  41),
                      cmap=cmaps.BlueWhiteOrangeRed,
                      transform=ccrs.PlateCarree())
    ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')

    ax4 = fig.add_subplot(3, 2, 4)
    ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
    ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')

    bb = ax4.bar(year, PCs[:, 1], color='r')
    for bar, height in zip(bb, PCs[:, 1]):
        if height < 0:
            bar.set(color='blue')

    ax5 = fig.add_subplot(3, 2, 5, projection=projection)
    mapart(ax5)
    ppp = ax5.contourf(file.lon,
                       file.lat,
                       EOFs[2],
                       np.linspace(-1,
                                   1,
                                   41),
                       cmap=cmaps.BlueWhiteOrangeRed,
                       transform=ccrs.PlateCarree())
    ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')

    ax6 = fig.add_subplot(3, 2, 6)
    ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
    ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')

    bbb = ax6.bar(year, PCs[:, 2], color='r')
    for bar, height in zip(bbb, PCs[:, 2]):
        if height < 0:
            bar.set(color='blue')

    # 为柱状图添加0标准线
    ax2.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax4.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax6.axhline(y=0,  linewidth=1, color='k', linestyle='-')

    # 在图下边留白边放colorbar
    fig.subplots_adjust(bottom=0.15)
    # colorbar位置： 左 下 宽 高
    l = 0.25
    b = 0.04
    w = 0.6
    h = 0.015
    # 对应 l,b,w,h；设置colorbar位置；
    rect = [l, b, w, h]
    cbar_ax = fig.add_axes(rect)

    c = plt.colorbar(
        pp,
        cax=cbar_ax,
        orientation='horizontal',
        aspect=20,
        pad=0.1) #绘制colorbar

    c.ax.tick_params(labelsize=14) #设置colorbar的标签大小

    plt.subplots_adjust(wspace=-0.1, hspace=0.3)
    plt.show()

eof_contourf(eof*(-1),pc*(-1),var,year)

D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)

north检验

north = solver.northTest(neigs=3).data
eig = solver.eigenvalues(neigs=3).data
print("north误差",north)
print("特征值",eig)

north误差 [115.06639   58.858727  25.045311]
特征值 [533.5409   272.9167   116.130325]

def northplt(north,eig,k):
    ispass = False
    for i in range (0,k-1):
        if (eig[i+1]+north[i+1])<(eig[i]-north[i]):
            print('第%d模态和第%d模态显著分离，通过north检验'%(i+1,i+2))
        else:
            print('第%d模态和第%d模态分离不显著，没有通过north检验'%(i+1,i+2))
            break
    else:
        ispass = True
    
    if ispass:
        plt.errorbar(np.arange(k)+1, eig,yerr=north)

        plt.xticks(np.arange(k)+1)
        plt.title('north test PASS')
        plt.show()

northplt(north,eig,3)

第1模态和第2模态显著分离，通过north检验
第2模态和第3模态显著分离，通过north检验

可以看到，三个主分量无重叠部分，显著分离

相关系数分析

from scipy.stats import pearsonr

def nonan(a):
    a = np.array(a)
    mask = np.isnan(a)
    a[mask] = 0
    b = np.reshape(a,(89*180))
    return b

def xiangguantest(eof,pc,k):
    for i in range(k-1):
        for j in range(i+1,k):
            r1,p=pearsonr(nonan(eof[i,:,:]),nonan(eof[j,:,:]))
            r1 = round(r1,2)
            print(f'第{i+1}模态和第{j+1}模态空间相关系数:{r1}')
            r2,p=pearsonr(pc[:,i],pc[:,i+1])
            r2 = round(r2,2)
            print(f'第{i+1}模态和第{j+1}模态时间相关系数:{r2}')
        
    

xiangguantest(eof,pc,3)

第1模态和第2模态空间相关系数:0.17
第1模态和第2模态时间相关系数:0.0
第1模态和第3模态空间相关系数:0.04
第1模态和第3模态时间相关系数:0.0
第2模态和第3模态空间相关系数:-0.02
第2模态和第3模态时间相关系数:-0.0

原始值

def mapart_(ax): #定义一个绘制地图地图的函数
    '''
    添加地图元素
    '''
    projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影，并指定中心经度为180°
    ax.coastlines(color='k', lw=0.5)                #添加海岸线
    ax.add_feature(cfeature.LAND, facecolor='white') #添加陆地
    # 设置地图范围，经度为（0，358），维度为（-88，88）
    ax.set_extent([0, 358, -88, 88], crs=ccrs.PlateCarree(central_longitude=180))
    # 设置经纬度标签
    ax.set_xticks([-120, -60, 0, 60, 120], crs=projection)
    ax.set_yticks([-80, -40, 0, 40, 80], crs=projection)
    lon_formatter = LongitudeFormatter(zero_direction_label=True) #给标签添加对应的N,S,E,W
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)


def eof_contourf_(EOFs, PCs, pers, year):  #定义一个绘制填充图的函数
    plt.close
    pers=(pers*100).values  #将方差转换为百分数的形式
    fig = plt.figure(figsize=(18, 12), dpi=300)       
    projection = ccrs.PlateCarree(central_longitude=180) 
    ax1 = fig.add_subplot(3, 2, 1, projection=projection)
    mapart(ax1) 
    p = ax1.contourf(file.lon,
                     file.lat,
                     EOFs[0], 

                     cmap=cmaps.BlueWhiteOrangeRed,
                     transform=ccrs.PlateCarree())  
    ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left') 

    ax2 = fig.add_subplot(3, 2, 2) #绘制第二个子图
    ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--') 
    b = ax2.bar(year, PCs[:, 0], color='r') 
    for bar, height in zip(b, PCs[:, 0]):
        if height < 0:
            bar.set(color='blue') 
    ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')

    ax3 = fig.add_subplot(3, 2, 3, projection=projection) 
    mapart(ax3)
    pp = ax3.contourf(file.lon,
                      file.lat,
                      EOFs[1],

                      cmap=cmaps.BlueWhiteOrangeRed,
                      transform=ccrs.PlateCarree())
    ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')

    ax4 = fig.add_subplot(3, 2, 4)
    ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
    ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')

    bb = ax4.bar(year, PCs[:, 1], color='r')
    for bar, height in zip(bb, PCs[:, 1]):
        if height < 0:
            bar.set(color='blue')

    ax5 = fig.add_subplot(3, 2, 5, projection=projection)
    mapart(ax5)
    ppp = ax5.contourf(file.lon,
                       file.lat,
                       EOFs[2],

                       cmap=cmaps.BlueWhiteOrangeRed,
                       transform=ccrs.PlateCarree())
    ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')

    ax6 = fig.add_subplot(3, 2, 6)
    ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
    ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')

    bbb = ax6.bar(year, PCs[:, 2], color='r')
    for bar, height in zip(bbb, PCs[:, 2]):
        if height < 0:
            bar.set(color='blue')

    # 为柱状图添加0标准线
    ax2.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax4.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax6.axhline(y=0,  linewidth=1, color='k', linestyle='-')

    # 在图下边留白边放colorbar
    fig.subplots_adjust(bottom=0.15)
    # colorbar位置： 左 下 宽 高
    l = 0.25
    b = 0.04
    w = 0.6
    h = 0.015
    # 对应 l,b,w,h；设置colorbar位置；
    rect = [l, b, w, h]
    cbar_ax = fig.add_axes(rect)

    c = plt.colorbar(
        pp,
        cax=cbar_ax,
        orientation='horizontal',
        aspect=20,
        pad=0.1) #绘制colorbar

    c.ax.tick_params(labelsize=14) #设置colorbar的标签大小

    plt.subplots_adjust(wspace=-0.1, hspace=0.3)
    # plt.tight_layout()
    # plt.savefig('./figure/eof2.png')
    plt.show()

# eof
solver = Eof(sst,center=False)
# solver = Eof(sst)
eof_ = solver.eofs(neofs=3,eofscaling=2)
pc_ = solver.pcs(npcs=3,pcscaling=1)
var_=solver.varianceFraction(neigs=3)

print('eof_:',eof_.shape)
print('-------------------')
print('pc_:',pc_.shape)

eof_: (3, 89, 180)
-------------------
pc_: (43, 3)

eof_contourf_(eof_*(-1),pc_*(-1),var_,year)

north检验

north_ = solver.northTest(neigs=3).data
eig_ = solver.eigenvalues(neigs=3).data
print("north误差",north_)
print("特征值",eig_)

north误差 [7.8232719e+05 8.7714035e+01 4.2884262e+01]
特征值 [3.6275018e+06 4.0671323e+02 1.9884615e+02]

northplt(north_,eig_,3)

第1模态和第2模态显著分离，通过north检验
第2模态和第3模态显著分离，通过north检验

相关系数分析

xiangguantest(eof_,pc_,3)

第1模态和第2模态空间相关系数:0.1
第1模态和第2模态时间相关系数:-0.59
第1模态和第3模态空间相关系数:0.12
第1模态和第3模态时间相关系数:-0.59
第2模态和第3模态空间相关系数:-0.02
第2模态和第3模态时间相关系数:-0.0

标准化

from sklearn  import  preprocessing
import eofs

按空间标准化

sst

sst_scaled_t = np.zeros((43,89,180))
for i in range(43):
    sst_scaled_t[i,:,:] = preprocessing.scale(sst.data[i,:,:])
    
sst_scaled_t

D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:235: UserWarning: Numerical issues were encountered when centering the data and might not be solved. Dataset may contain too large values. You may need to prescale your features.
  warnings.warn(
D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:254: UserWarning: Numerical issues were encountered when scaling the data and might not be solved. The standard deviation of the data is probably very close to 0. 
  warnings.warn(

array([[[-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
         -1.29350352, -1.29724193],
        [-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
         -1.29350352, -1.29724193],
        [-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
         -1.29350352, -1.29724193],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]],

       [[-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
         -1.30991757, -1.31138611],
        [-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
         -1.30991757, -1.31138611],
        [-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
         -1.30991757, -1.31138611],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]],

       [[-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
         -1.29943419, -1.30155838],
        [-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
         -1.29943419, -1.30155838],
        [-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
         -1.29943419, -1.30155838],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]],

       ...,

       [[-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
         -1.33291507, -1.33418429],
        [-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
         -1.33291507, -1.33418429],
        [-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
         -1.33291507, -1.33418429],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]],

       [[-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
         -1.32402062, -1.32561731],
        [-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
         -1.32402062, -1.32561731],
        [-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
         -1.32402062, -1.32561731],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]],

       [[-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
         -1.3279072 , -1.32902801],
        [-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
         -1.3279072 , -1.32902801],
        [-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
         -1.3279072 , -1.32902801],
        ...,
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan],
        [        nan,         nan,         nan, ...,         nan,
                 nan,         nan]]])

solver = eofs.standard.Eof(sst_scaled_t)
eof__ = solver.eofs(neofs=3,eofscaling=2)
pc__ = solver.pcs(npcs=3,pcscaling=1)
var__=solver.varianceFraction(neigs=3)
print('eof__:',eof__.shape)
print('-------------------')
print('pc__:',pc__.shape)

eof__: (3, 89, 180)
-------------------
pc__: (43, 3)

eof_contourf_(eof__*(-1),pc__*(-1),var_,year)

D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)

north检验

north__ = solver.northTest(neigs=3)
eig__ = solver.eigenvalues(neigs=3)
print("north误差",north__)
print("特征值",eig__)

north误差 [0.51234382 0.24994451 0.15094646]
特征值 [2.37564056 1.15894504 0.69990995]

northplt(north__,eig__,3)

第1模态和第2模态显著分离，通过north检验
第2模态和第3模态显著分离，通过north检验

相关系数分析

xiangguantest(eof__,pc__,3)

第1模态和第2模态空间相关系数:-0.0
第1模态和第2模态时间相关系数:0.0
第1模态和第3模态空间相关系数:0.0
第1模态和第3模态时间相关系数:0.0
第2模态和第3模态空间相关系数:0.0
第2模态和第3模态时间相关系数:0.0

时间序列标准化

sst

<xarray.DataArray 'sst' (time: 43, lat: 89, lon: 180)>
array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        ...,
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]],

       [[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
...
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]],

       [[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        [-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
         -1.7999998, -1.7999998],
        ...,
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
  * lon      (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
  * time     (time) int64 1979 1980 1981 1982 1983 ... 2017 2018 2019 2020 2021

sst_scaled_p = np.zeros((43,89,180))
for i in range(89):
    for j in range(180):
        sst_scaled_p[:,i,j] = preprocessing.scale(sst.data[:,i,j])

sst_scaled_p

D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:235: UserWarning: Numerical issues were encountered when centering the data and might not be solved. Dataset may contain too large values. You may need to prescale your features.
  warnings.warn(
IOPub data rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_data_rate_limit`.

Current values:
ServerApp.iopub_data_rate_limit=1000000.0 (bytes/sec)
ServerApp.rate_limit_window=3.0 (secs)

solver = eofs.standard.Eof(sst_scaled_p)
eof__p = solver.eofs(neofs=3,eofscaling=2)
pc__p = solver.pcs(npcs=3,pcscaling=1)
var__p=solver.varianceFraction(neigs=3)
print('eof__p:',eof__p.shape)
print('-------------------')
print('pc__p:',pc__p.shape)

eof__p: (3, 89, 180)
-------------------
pc__p: (43, 3)

def eof_contourf__(EOFs, PCs, pers, year):  #定义一个绘制填充图的函数
    plt.close
    pers=(pers*100)  #将方差转换为百分数的形式
    fig = plt.figure(figsize=(18, 12), dpi=300)       
    projection = ccrs.PlateCarree(central_longitude=180) 
    ax1 = fig.add_subplot(3, 2, 1, projection=projection)
    mapart(ax1) 
    p = ax1.contourf(file.lon,
                     file.lat,
                     EOFs[0], 

                     cmap=cmaps.BlueWhiteOrangeRed,
                     transform=ccrs.PlateCarree())  
    ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left') 

    ax2 = fig.add_subplot(3, 2, 2) #绘制第二个子图
    ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--') 
    b = ax2.bar(year, PCs[:, 0], color='r') 
    for bar, height in zip(b, PCs[:, 0]):
        if height < 0:
            bar.set(color='blue') 
    ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')

    ax3 = fig.add_subplot(3, 2, 3, projection=projection) 
    mapart(ax3)
    pp = ax3.contourf(file.lon,
                      file.lat,
                      EOFs[1],

                      cmap=cmaps.BlueWhiteOrangeRed,
                      transform=ccrs.PlateCarree())
    ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')

    ax4 = fig.add_subplot(3, 2, 4)
    ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
    ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')

    bb = ax4.bar(year, PCs[:, 1], color='r')
    for bar, height in zip(bb, PCs[:, 1]):
        if height < 0:
            bar.set(color='blue')

    ax5 = fig.add_subplot(3, 2, 5, projection=projection)
    mapart(ax5)
    ppp = ax5.contourf(file.lon,
                       file.lat,
                       EOFs[2],

                       cmap=cmaps.BlueWhiteOrangeRed,
                       transform=ccrs.PlateCarree())
    ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')

    ax6 = fig.add_subplot(3, 2, 6)
    ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
    ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')

    bbb = ax6.bar(year, PCs[:, 2], color='r')
    for bar, height in zip(bbb, PCs[:, 2]):
        if height < 0:
            bar.set(color='blue')

    # 为柱状图添加0标准线
    ax2.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax4.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax6.axhline(y=0,  linewidth=1, color='k', linestyle='-')

    # 在图下边留白边放colorbar
    fig.subplots_adjust(bottom=0.15)
    # colorbar位置： 左 下 宽 高
    l = 0.25
    b = 0.04
    w = 0.6
    h = 0.015
    # 对应 l,b,w,h；设置colorbar位置；
    rect = [l, b, w, h]
    cbar_ax = fig.add_axes(rect)

    c = plt.colorbar(
        pp,
        cax=cbar_ax,
        orientation='horizontal',
        aspect=20,
        pad=0.1) #绘制colorbar

    c.ax.tick_params(labelsize=14) #设置colorbar的标签大小

    plt.subplots_adjust(wspace=-0.1, hspace=0.3)
    # plt.tight_layout()
    # plt.savefig('./figure/eof2.png')
    plt.show()

eof_contourf__(eof__p*(-1),pc__p*(-1),var__p,year)

D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)

north检验

north__p = solver.northTest(neigs=3)
eig__p = solver.eigenvalues(neigs=3)
print("north误差",north__p)
print("特征值",eig__p)

north误差 [686.41866511 301.32724554 152.68502965]
特征值 [3182.79241419 1397.1969587   707.97135746]

northplt(north__,eig__,3)

第1模态和第2模态显著分离，通过north检验
第2模态和第3模态显著分离，通过north检验

相关系数分析

xiangguantest(eof__p,pc__p,3)

第1模态和第2模态空间相关系数:0.2
第1模态和第2模态时间相关系数:0.0
第1模态和第3模态空间相关系数:-0.01
第1模态和第3模态时间相关系数:0.0
第2模态和第3模态空间相关系数:0.01
第2模态和第3模态时间相关系数:-0.0